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Abstract 

The simulation of the metabolism in mammalian cells becomes a severe problem if spa- 
^ ■ tial distributions must be taken into account. Especially the cytoplsma has a very complex 

geometric structure which cannot be handled by standard discretization techniques. In the 
present paper we propose a homogenization technique for computing effective diffusion con- 
stants. This is accomplished by using a two-step strategy. The first step consists of an ana- 
£Sj . lytic homogenization from the smallest to an intermediate scale. The homogenization error 

is estimated by comparing the analytic diffusion constant with a numerical estimate obtained 
by using real cell geometries. The second step consists of a random homogenization. Since 
no analytical solution is known to this homogenization problem, a numerical approxima- 
tion algorithm is proposed. Although rather expensive this algorithm provides a reasonable 
^ ' estimate of the homogenized diffusion constant. 
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1 Introduction 

When mammalian cells are exposed to foreign and potentially harmful compounds a series of 
events takes place. Following uptake the substance is distributed in different intracellular com- 
partments by diffusion, absorption and desorption. The majority of the compound is either dis- 
solved in the aqueous phase, the cytoplasm, or in the lipophilic phase, the membranes. Parallel 
to diffusion and absorption/desorption bioactivation/biotransformation by different soluble and 
membrane bound enzymes takes place. The purpose of biotransformation is to render the sub- 
stance suitable for excretion. 

A human cell consists schematically of an outer cellular membrane, a cytoplasm containing 
a large number of organelles (mitochondria, endoplasmatic reticulum etc.), a nuclear membrane 
and finally the cellular nucleus containing DNA. Figure Q] shows a sketch of a cell while Figure[2] 
shows a microphotograph of a nucleus with part of the surrounding cytoplasm. The organelle 
membranes create a complex and dense system of membranes or subdomains throughout the 
cytoplasm. The mathematical description of the biotransformation leads to a system of reaction- 
diffusion equations in a complex geometrical domain, dominated by thin membranous structures 
with similar physical and chemical properties. If these structures are treated as separate sub- 
domains, any model becomes computationally very expensive. Moreover, due to the natural 
variations in the cell structures, every individual cell needs its own mathematical model. 

In order to make the system numerically treatable while at the same time retaining the essen- 
tial features of the metabolism under consideration, in 01 a way of homogenizing the cytoplasm 
has been developed, aiming at a manageable system of reaction-diffusion equations for the vari- 
ous species. In the present paper, we report about numerical experiments which justify some of 
the strategies in the cited paper. The general modelling assumptions are summarized below. We 
will use them also in the present report. 

Modelling assumptions: 

• On a small scale in space, the volume between the outer cellular membrane and the nucleus 
membrane consists of layered structures cytoplasm/membranes. 

• In the large scale, this volume contains an unordered set of the small-scale substructures 
which are uniformly distributed over the volume. 



3 



Cell Structure 
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Figure 1: Schematic picture of a cell. Picture copyright U.S. National 
Cancer Institute's Surveillance, Epidemiology and End Results Program, 
http : //training, seer . cancer .gov/module_anatomy/unit2_l_cell_functions_l .html 




Figure 2: Ultrastructure of the cell, nucleus and cytoplasm. Picture copyright Histology Learning 
System, Boston University, http : //www . bu . edu/histology/m/index . htm 
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2 THE MATHEMATICAL MODEL 



• The physical and chemical properties of the cytoplasm and of the membranes are uniform. 

• We adopt the continuum hypothesis, i.e., we assume that the set of molecules in the cell 
can be modelled by considering a continuous representation (a concentration). 

• The processes of absorption and desorption of the individual species into or out of the 
membrane is much faster than the diffusion and reaction processes. In this case, the re- 
lations of the concentrations of a species near a membrane/cytoplasm boundary can be 
conveniently described with the help of a partition coefficient. 

In Section [2l we introduce the mathematical model. The next section is devoted to a descrip- 
tion of the general experimental set-up which is used to compute effective diffusion coefficients 
numerically. For the solution of the arising boundary-value problems for partial differential 
equations we used the Comsol Multiphysics^l [2] environment. 

In [?], the diffusion coefficients in the membrane structures have been homogenized by as- 
suming the membranes and the aqueous volumes to be ideal infinite plane layers. This allows 
for an analytic computation of the effective diffusivity. In Section SI we will compare effective 
diffusivities obtained this way against numerically determined effective diffusivities by using 
computational domains which have been discretized from microphotographs of cell membranes. 

The result of the first homogenization step leads to anisotropic diffusion tensors valid locally. 
Invoking the assumption about the random distribution of the orientation of the membranes, the 
next step consists of a stochastic homogenization. In contrast to the one- and two-dimensional 
case, no analytical solution in the general three-dimensional case is known. We will compute the 
effective diffusion coefficient numerically by Monte Carlo techniques in Section [51 

2 The Mathematical Model 

2.1 The Governing Equations 

We intend to derive a homogenized model of the reaction and diffusion processes inside the 
cytoplasm. For that purpose, let G denote the volume between the outer cell membrane and the 
nucleus membrane (excluding the membranes themselves). This volume is split into two disjoint 
parts G; and G w which denote the lipophilic part and the aqueous part, respectively, of the cell. 
Note that these subdomains are not necessarily connected. Assume that we are interested in the 
contrations ci, . . . ,c n of n species inside of G. For the k-th species, it holds 

—c k = V-(d k (x)c k )+R k (ci,...,c n ,x), xeG, k=l,...,n. (1) 

Here, d k denotes the diffusion tensor of the k-th species which is assumed to be constant in both 
G; and G w . R k denotes the reaction term. It varies strongly with x. In the lipophilic part, R k = 
because no reactions are taking place there. The concentrations of some of the species can be 
assumed to be constant over time. As a consequence, many of the reaction terms will be linear. 



Comsol Multiphysics is a registered trade mark of Comsol AB, Stockholm, Sweden. 



2.2 The Model Problem 
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The partition coefficient, K pk is the equilibrium ratio of the concentrations of species k be- 
tween the aqueous compartment and its adjacent lipid compartment. This gives rise to boundary 
conditions, 

Ck,w = K Pik c kj i * G G w fl G/, (2) 

on the inner boundaries where c^ w and c^i denote the concentrations in the aqueous and lipid 
parts, respectively. 

The system (OQ) with inner boundary conditions © will be supplemented by (outer) boundary 
conditions and initial conditions. For the purposes of this paper the precise structure of these 
conditions is not important. 



2.2 The Model Problem 



Let G C R 3 be a bounded domain which will be splitted into two (not necessarily connected) 
subdomains G\ and G2 such that 

GiHG 2 = 0, G = GiUG 2 . 
The interior boundary will be denoted by T, 

T = GiHG 2 . 

Consider the equations 
d 



dt 



Vj - V ■ (di(x)Vvj) + n(x)vi = fi(x) , xeGi, i = l, 2. 



(3) 



Assume additionally boundary conditions on dG and initial conditions on G be given. 

On the inner boundary T the flux must be continuous. Let n ; - denote the outer normal at the 
boundary of G ; . The continuity conditions reads now 



\-U2- 



0, xeT. 



(4) 



dni ' " z dii2 

The presence of a partion coefficient between the two phases gives rise to the boundary condition 

v\=K p V2, xeT. (5) 
This problem can be reduced to a problem in a more standard form by introducing 



u(x) 



fvi(x), xeG\, 
[K p v 2 (x), xeG 2 . 



(6) 



For this new function u, the inner boundary conditions become 

du 



<9ni 



1 , du 
+ -^d 2 



dn 2 



0, 



G 2 



xeT. 



lev 
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2 THE MATHEMATICAL MODEL 



This motivates the definitions 

d _\di, xgGi, a _f 1 ' xeGu 
\d 2 /K p , xeG 2 , \l/Kp, xeG 2 , 

r _j ri ' xGG i' /=|^ 1 ' xGGl ' 

With these definitions, the problem ©, © becomes equivalent to 

d 

a^-u-V ■(dVu) + ru = f, x<EG (7) 

subject to correspondingly modified initial and boundary conditions. 
For later use, let 

\ G l\ i \ G 2\ 



2.3 Going From The Smallest to The Medium Scale: Homogenization of a 
Periodic Structure 

The homogenization procedure for an equation of the type (Q is proved in [6]. We cite the basic 
facts. Consider the following problem, 

o e -u e +A e u s = f, u £ (0) = u 

dt (9) 

u e eL 2 (0,T;H^{G)). 

Here, the operator A 8 is given by 

A £ u := ( dfj^—u] + r £ u. 

OXi \ J OXj J 

For convenience, we use the Einstein summation convention: If an index appears twice in a 
multiplicative expression, this expression is understood to implicitly represent the sum over this 
expression where the index varies between 1 and 3 (the dimension of G). Moreover, we assume 
the following construction of the coefficients: 

o £ {x) = a(x/e), r £ {x) = r{x/e), dfj(x) = dij(x/e), ij =1,2,3. 

The functions o £ and r £ are assumed to belong to L°°(G), and 

o > Co > 0, r{x) > a.e. in G 

for some Oq E R. The functions dfj are assumed to be measurable and to satisfy the conditions 
dfi = d% and 

a\E, | 2 < dfj&£j </3|^| 2 , a.e. in G for all £, e R 3 and < a < /3 < oo. 



2.3 Going From The Smallest to The Medium Scale: Homogenization of a Periodic Structure! 



Finally, assume u Q e L 2 (G) and f e L 2 (0, T;L 2 (G)). 
In order to find the homogenized equation, assume 

f — weakly in L 2 (0, T;L 2 (G)). 

Let Y be an axis parallel hexahedron in M 3 , that is, 

3 

Y = x (a;,fcj). 
i=i 

For a y-periodic function /, the mean value is given by 

(/) := ^ / 'f{y)dy. 

Y 

Assume now that a*/, <7, and r are y-periodic. Then it is possible to consider the problem 

(o)—u+Au = f, u(0) = u 
ueL 2 (0,T;H^(G)), 



where the operator A is given by 



d ( . d 



Au = -— ld eS jj—u ) + (r}u 



I d( Pj\ 
«eff,(7 = (dij — "ik-^—^ ) j 

and <Pj is the y-periodic solution of the following local elliptic problem: 



(PyeW(y). 



Here, W(7) = {0 G J ff 1 (y)|(pis y-periodic and (<p) = 0}. 

Theorem 2.1. Under the conditions stated above, and rfTOl) /zave unique solutions u E G 
L 2 (0, r;i/Q (G)) and m 6 L 2 (0, ^;#o (^))> respectively, and it holds 

u £ — > u in L 2 (0, T;Hq (G)) weakly as £ — > 0. 

This theorem is proved in [6, p. 56] H 

In our model problem © the cell problems can be simplified considerably. We will assume 
that, in the smallest scale, aqueous and lipid compartments are perfectly layered. It turns out 



In the reference, the proof is given for a problem without reaction term. But the proof can easily be generalized. 
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2 THE MATHEMATICAL MODEL 



that in this case the computation for the effective diffusion coefficients leads to a transmission 
problem which can be solved analytically. 

Consider the cell problem (fTTj) . The material consists of perfect layers co + and CO with 
thicknesses a + and a~, respectively. The diffusion coefficients in these layers are d + and d~, 
respectively. According to our assumptions, 

— - El 

Let us choose the following coordinate: x 2 is normal to the interface between co + and CO while 
x\ and xt, span this interface. The diffusion coefficient d is given by 





[0, 


if / ^ j, 


dij(x) = | 




if i = j and x G CO + 




W, 


if i = j and x G C0~ 



The cell problem is posed on 

Y = (0 i h)x{-a-,a + )x(0,h). 
Then the homogenized diffusivities become (see Q) 

^eff,n = (a + d^ +a~d^)/(a + + a~), (12) 

<4ff,22 = ( a+ + a ~)/ { a+ 1&2 + a /d 2 ) i (13) 

4ff,33 = [a + d% +a'd^)/(a + +a~), (14) 

dij = if i^j. (15) 

Note that J e ff,n an d <i e ff,33 are the arithmetic means while rf e ff,22 is the harmonic mean of both 
diffusivities a\ and aj . 

2.4 From The Medium to The Large Scale: Stochastic Homogenization 

In global coordinates, we cannot assume that the coordinate system is oriented in the way that 
we used above. Consider two Cartesian coordinate systems {x\,X2,xj) and (zi,Z2>£3)- Assume 
that a given point x has the representation z = Tx with respect to the z-coordinates. Note that T 
is an orthogonal matrix in that case. Denote the matrix of diffusion coefficients with respect to 
the x-coordinates by Q* and that with respect to the z-coordinates by Q. Then a short calculation 
yields 

Q* = TQT 1 = TQT T . 

This is the point to invoke the next critical assumption: We assume that the volume is tightly 
packed with substructures of the type considered before, namely layered materials. The key 
assumption is that all orientations are equally probable. This leads to a stochastic description of 



2.4 From The Medium to The Large Scale: Stochastic Homogenization 
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the diffusion coefficients. We need a homogenization of operators with random coefficients. A 
theory for that is provided in |0. 

Note that the mean values (a) and (r) in (flOl ) are independent of the orientation of the layers. 
Therefore, it suffices to consider the stationary diffusion problem 

A £ u £ =f 1 u £ EH^(G) (16) 

with G c R m and 

Aeu= -^{ dfj i u ) 

which is the counterpart of ©. Assume as before that 

d fj ( x ) = d ij ( x / £ ) > h j = 1 , • ■ ■ , m. 

The randomness of the orientation is modelled by assuming that the matrix A(y) = (dij(y)) is 
statistically stationary with respect to the spatial variable y E W n , or equivalently, that A(y) is a 
typical realization of a stationary random field. 

Let (Q,,J?,P) be a probability space with c-algebra & and probability measure P. Let for 
each x 6 W" a random variable £ (x) over (Q, ^,P) be given. The random field t, is stationary if 
it can be represented in the form 

£(x, co) = a{T{x)(o) 

where a(-) is a fixed random variable, T = T(x) : Q — > SI is a measurable transformation which 
preserves the measure P on (£2, J^) . Therefore, for the definition of the coefficients dtj in (fT6l) it is 
sufficient to consider a matrix ( ) of random variables : H — > R. Realizations of coefficients 
can then be obtained by setting 

dij(x,(o) =dij(T(x)(D). 
Assume in the following that J 7 y G L°°(Q) and 

a|^| 2 <J ;7 (co)^, ^Gl m 

for almost all co e Q. with a > independent of £ and ft). 

A (deterministic) matrix d(y) is said to admit a homogenization if there exists a constant 
elliptic matrix J e ff such that for any / G H l (G) the solutions w e of the Dirichlet problem (fT6l) it 
holds 

w e — )■ w in //J (G) weakly as £ — 0, and 
d £ Vu £ — > d s ffVu in L 2 (G) weakly as £ — > 0, 

where u is the solution of the Dirichlet problem 

-V-(deffV«)=/, ^^(G). 

This definition correspondents to the stationary version of Theorem |2.11 The following theorem 
holds true @p. 230]: 



10 



3 NUMERICAL DETERMINATION OF EFFECTIVE DIFFUSIVITIES 



Theorem 2.2. Assume additionally that the family of mappings T(x) : £1 — > £1, x G W n , forms an 
ergodic m-dimensional dynamical system. Then for almost all CO&Q., the matrix with coefficients 
dij(x) = dij(T(x)(o) admits homogenization, and the homogenized matrix d e jf is independent of 
CO. 

Unfortunately, an analytical representation of d e ff is only possible in exceptional cases. We 
are interested in diffusion coefficients having a representation 

d[x) = T(x,a)QT(x,G))~ 1 

where T(x) E SO(m) is uniformly distributed in SO(m) and Q is a fixed diffusion tensor. In the 
two-dimensional case, an analytic solution is provided in [|51 p. 235]. For m = 2, d e ff is simply a 
scalar equal to the geometric mean of the eigenvalues of Q, 

Here det(<2) denotes the determinant of Q. 

There is no analytical solution known for the case m = 3. 

For later use in the experimental estimation of the effective diffusivity, the following obser- 
vation is important: According to our assumptions on d, the estimate 

II (AT 1 ! < cr 1 

holds true such that, for any / and / in H~ 1 (G), 

\(l,(AT l f)\<\\l\\H-HG)U- l \\f\\H-i(G) 

independently of CO € £1. Consequently, for the expectation values it holds 

E(l,(A £ )- l f)^E(l,A- l f) (17) 
by the dominated convergence theorem. 

3 Numerical Determination of Effective Diffusivities 

Under the assumption that an effective diffusivity for a given problem exists, the corresponding 
diffusion constants can be determined experimentally. For that, let D C G be a subdomain which 
is in size comparable to G such that the small scale structure is considerable smaller than the size 
of D. Assume that we want to determine the (scalar) diffusion constant for the diffusion process 
in x-direction. In that case it is convenient to use a cylindrical domain 

D= (0,L) x co 

with ft) C M 2 being some bounded domain. On D consider the stationary diffusion equation 

-V • (d(x)Vu) = 0, xeD. 
The boundary conditions are selected as follows: 
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• On the boundary Tq = {0} x 0), a fixed Dirichlet condition is given, 

"r = c o- 

• On the boundary T\ = {1} x (0, a free diffusion into the surrounding medium is assumed, 

— n ■ (d(x)Vu)Y l = M(«n — c\). 

Here, M is the mass transfer coefficient and c\ is the concentration in the bulk solution 
outside of D. 

• All other boundaries 1^2 = dG\ (Tq UT\) are isolated, 

— n ■ (d(x)Vu)r-j = 0. 
If d(x) would be a constant J e ff, it would hold 

d rr c °~ Mout - N 

"erf — iv average 5 



^average = jjTT J M(u-Ci)dF, 

1 r. 



"out = jp-r J udT. 

By \Y\ I we denote the Lebesgue measure of T\. If d(x) is varying, these equations can be used 
as an estimation of the effective diffusivity d e ff . In case of an anisotropic effective diffusivity, the 
above construction leads to an estimate of the effective diffusivity in x-direction, i.e., J e ffji- 

In the one-dimensional case m = 1, however, an analytic solution is possible. A simple 
calculation gives 



L * - 1 



d e ff z 

which amounts to the harmonic mean. 




4 Theoretical And Experimental Diffusivities For Layered Struc- 
tures 

The homogenization of layered structures in Section 12.31 made use of the assumption that we 
have ideal planes of different materials with different diffusion tensors. In a real biological cell, 
this assumption is only approximately fulfilled in small subdomains. Besides the effect of not 
having the parameter e close to zero an additional error is introduced this way. The aim of the 
present section is to obtain some experimental estimates of how large the error will be. We will 
start with a real photograph of some cell organelles and extract the geometrical structure of the 
lipophilic and aqueous layers. Then the diffusivity is estimated using the strategy of Section [3j 
This diffusion constant will be compared to the theoretical homogenized value. 
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Figure 3: Detail of a rat cell showing the Golgi-apparatus. The box indicates the area used as a 
reference domain. Copyright Dr. H. Jastrow 

4.1 The Experimental Set-up 

The experiments of this section are based on the micro-photograph shown in Figure [3j The 
part enclosed by a box in that figure has been extracted and amplified in contrast. This way, the 
membrane structure in Figure|4]has been obtained. Note that only the black lines represent mem- 
branes. The geometry of this structure has been too complex for the software used in the numer- 
ical experiments. The number of degrees of freedoms obtained after discretization has become 
too large. Therefore, we extracted again a part of this geometry in order to make the problem 
tractable with the available software. Note that the diffusion in this problem is anisotropic. In 
order to be able to compare the experimental numerical diffusivity with the analytical value, the 
main orientation of the membranes was aligned with the y-axis. The resulting geometries can be 
found in Figure [51 Two cases have been considered. 

• Case A: In this case, almost perfect layers have been used. 

• Case B: Here we want to estimate the influence of short circuits and more irregular struc- 
tures. 

The geometrical data for both data are provided in Table [Q The corresponding data for 
the diffusion constants are given in Table [21 Observe that the diffusion in the lipophilic part is 
anisotropic. This has been used for the numerical experiments. In contrast to that, the homoge- 
nized diffusion constant has been determined by using d2,n> only. So we expect a larger error in 
the experiments with the domain of case B. 



The Experimental Set-up 




Figure 4: Contrast amplified reference domain. The black areas indicate membranes 




























(a) 



I Mil 

(b) 

Figure 5: Computational domains for case A (a) and case B (b) 





Case A 


CaseB 




4.359 x 1CT 7 


4.390 x 10" 7 


l y [m] 


0.9125 x 1(T 7 


1.568 x 10- 7 


Pi 


0.8122 


0.8139 


P2 


0.1878 


0.1861 



Table 1 : Geometric constants 
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Value 


di[m z s- y ] 
^.llfmV 1 ] 
J 2 ,22[m 2 s _1 ] 
K P 


1.0 x 10" i4 
1.0 x 10~ 12 
1.0 x 10" 10 
1.26 x 10" 2 



Table 2: Diffusion constants 



The experimental determination of the effective diffusivity according to Section [3] can be 
carried out using two approaches: 

1. Use the original equation © subject to the inner boundary conditions @i and ©. 

2. Use the transformed problem © without any inner boundary conditions. 

In order to be as close as possible to the original problem we have chosen the first alternative for 
our experiments. Note that, in the case of K p = 1, both approaches are identical. 

Unfortunately, it is not possible to formulate the inner transmission conditions ©, © directly 
in Comsol Multiphysics. Instead, both conditions have been coupled by a penalty approach as 
suggested in [3]. For a suitably chosen constant K, ©, ©) is replaced by 

dv 

di^- = K(vi-K p V2), xeT (inGi), 

^ (18) 
d 2 -^ = k(K p v 2 -vi), xeT (inG 2 ). 

K acts as a mass transfer coefficient. 

For comparison purposes, even the homogenized problem (fTOT ) has been implemented in 
Comsol Multiphysics. 



4.2 Results 

The experiments have been carried out using the values 

k=10~ 4 , M=10~ 7 , c = l, ci=0. 

The penalty parameter has been chosen such that both sides of the equations (TT8T) are somehow 
in balance. The value of M has been chosen such that the outflow has the order magnitude 0.3co- 
In case 1, K p = 1 while, in case 2, K p = 0.0126. The results are summarized in Tabled 

The effective diffusivities given above refer to the steady state. In order to get a feeling for 
the influence of the homogenization on the transient behavior, we compared the time history of 
the mean flux out of the domain at the left boundary between the original equation ([3]) and its 



4.2 Results 
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case 


horn, constant 


exper. constant 


rel. difference 


1A 


1.2284 


1.3131 


6.9% 


IB 


1.2258 


1.3590 


10.9% 


2A 


1.2312 


1.2910 


4.9% 


2B 


1.2286 


1.4680 


19.5% 



Table 3: Homogenized and experimental effective diffusivities scaled by 10 14 
x i o -8 Avarage flux at the outflow boundary 



3.5 



3 - 




Figure 6: Comparison of the flux at the outflow boundary for the homogenized model (line) and 
the detailed model (dashed line) 

homogenized counterpart (flOl) . For that, the boundary value problem has been solved as before 
using the initial condition 

u(x,y) = l(T 6 exp ^-1000 ^ 1Q _ ? ) j aU = 0. 

The value of the experimental effective diffusion has been used in the homogenized problem. 
The results for the four different cases differ only marginally. As expected from the experiments, 
the largest differences occur in case 2B. This is shown in Figure [6] 
Summarizing, the following sources of errors occur: 
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5 EXPERIMENTAL EFFECTIVE DIFFUSIVITIES IN RANDOM MEDIA 



• The sizes of the sub-structures are not infinitesimal small; 

• The membrane layers are not ideal planes; 

• In the geometry case B, the membranes are touching the outflow boundary; 

• For computing the homogenized diffusion coefficient, only the normal part of the mem- 
brane diffusion tensor has been used; 

• The partition coefficients are handled by a penalty approach. 

The size of the numerical errors is negligible compared to the ones given above. 

5 Experimental Effective Diffusivities in Random Media 
5.1 The Experimental Set-up 

The idea for estimating the effective diffusivity in the present case is to use a Monte Carlo 
simulation. For that, the test domain D of Section |3] is chosen to be the unit cube, D = (0, l) 3 . 
Let 

Q = diag(Ji 1,^22,^33) 
be a fixed diffusion tensor. For a given N e N, this cube is subdivided into N 3 sub-cubes 

Dijk = (xi-i,Xi) x (yj-uyj) x (z k -uz k ) 

with Xj = y, = n = ih and h = 1 /N. h plays the role of £ in Theorem 12 .21 One experiment consists 
of choosing a realization d £ such that 

d£ \D ijk = T ijkQTij k 

where Ty* £ 50(3) are drawn uniformly distributed in 50(3). 

In order to describe the orientation we will use the Euler angles. Any rotation in 50(3) can 
be described by three angles, the so-called Euler angles. We will use the convention to first rotate 
around the *3-axis by the angle a, then around the (new) jq-axis by j8, and finally around the 
new ^3-axis by 7. This can be described formally by 

T=R 3 {y)R l {P)R 3 {a), a, ye (0,2*), j3<E(0,rc), (19) 

where 

(cos \j/ sin \j/ 0\ / 1 

-sinxj/ cos xjf , /?i(/3) = cos/3 sin/3 

1/ \0 -sin/3 cos/3. 

Let /x denote the Haar measure on 50(3). Its density has the simple form 

d\L = sin 3d ad fidy 
8tt z 



5.2 Results in 2D 
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with respect to the Lebesgue measure on (0, 2n) x (0, n) x (0, 2n) . 

This way, the expectation value of J| ff can be estimated for given /V (£). The computation of 
Wout and Average consists essentially of the evaluation of integrals 

J udT 
r, 

for functions u E H l (G). Since this is a continuous linear functional, we obtain 

Jf ff — > Jeff for e — > 

by using (fTTT) . Since there are no preferred directions in this setting, the diffusion is isotropic 
such that J e ff is a scalar. 

5.2 Results in 2D 

In the two-dimensional setting, an analytic solution of the random homogenization problem is 
known. Let 

Q = diag(J 1 i, d 2 2). 
The effective diffusivity is the scalar (51 p. 235] 

4ff = (JllJ22) 1/2 - 

We will carry out the experiment described above in the two-dimensional setting in order to 
obtain a certain gauge for its three-dimensional equivalent. 

The two-dimensional counterpart of the experiment described in Section 15.11 is to choose 
D = (0, l) 2 which will be subdivided, for a given N EN, into sub-squares 

Dij = (xi-uxj) x (yj-uyj) 

with Xi = yt = ih and h = l/N. The realizations d £ are now described by 

d e \ Dij = TijQT^ 

where 7}j E 50(2) are sampled uniformly distributed in 50(2). The elements of 50(2) are simple 
rotations described uniquely by an angle <p E [0, 2k), 

_ ( cos<p sin<p\ 
sin<p cos (p J ' 

The Haar measure ji on 50(2) has the density d\l = j^d(p with respect to the Lebesgue measure 
on (0,27t). 

The experimental results for 

Jn = 1, J22 = 10, J e ff = 3.1623 

are provided in Table HI We can draw the following conclusions: 
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6 CONCLUSIONS 



• The main parameter for the accuracy of the estimation of the effective diffusivity is N. This 
isn't hardly surprising. 

• For a given N, the sample size has only a minor influence on the accuracy. Once a certain 
number of trials has been reached the accuracy does not become better. The optimal sample 
size seems to be independent of N. 

• The standard deviation for sufficiently large sample sizes roughly halves while doubling 
N. This indicates a linear rate of convergence. 

• In all experiments, the mean value of the experimental effective diffusivity is an overesti- 
mation of the exact value. 

• If the sample size is too small, the standard deviation is misleading small. 

• The experiments in Section |4] indicated an error in the order of magnitude of 5% - 10% 
between the theoretical homogenized diffusivity and the experimentally observed. These 
results suggests to use a value of N = 20 and a sample size of at least 15 trials. 



5.3 Results in 3D 

Finally, the experimental estimation in the three-dimensional case has been carried out. Unfor- 
tunately, the geometry handling in Comsol Multiphysics has led to a severe restriction on how 
large N can be. Although the machine used had enough memory installed (16 GB), the Java 
heap space got exhausted rather soon. Moreover, the geometry analysis was surprisingly time- 
consuming compared to the assembly and solution process. 

For the experiments, we have chosen the diffusion coefficients 

Jil =9, d 2 2 = 25, J 33 = 1. 

Note that an analytical solution of the homogenization problem is not known. The results are 
given in Table [51 

6 Conclusions 

The present paper explains the homogenization strategy which has been used to derive effective 
equations for modelling the detailed metabolism in mammalian cells. The cytoplasm has been 
modelled assuming that three different length scales can be observed. For going from the smallest 
to the medium scale, an analytic homogenization technique is used. By comparing the analytic 
effective diffusion constant with results from numerical simulations on real cell geometries taken 
from photographs an error of 5% - 20% has been observed. Given the accuracy of the known 
diffusion constants in the lipophilic and aequous parts of the cytoplasm this accuracy appears to 
be sufficient. 
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Table 4: Experimental effective diffusivities in 2D for Jn = 1, J22 = 10, J e ff = 3. 1623 
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Table 5: Experimental effective diffusivities in 3D for d\\ =9, J22 = 25, J33 = 1 
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For the step from the medium scale to the large scale, a random homogenization technique 
has been used. Matheatically the effective diffusivity is known to exists. In the present paper an 
algorithm has been developed and tested for estimating the homogenized diffusion constant on 
the large scale. However, the computation times in Comsol Multiphysics have become very large 
(up to one week on a compute server based on a 2GHz AMD Opteron processor) for a reasonable 
setup such that alternative solution techniques should be investigated. 

The critical assumption in the last step is that about the probability distribution of the struc- 
tures on the intermediate scale. Its validity can probably only be justified by comparision to 
biochemical experiments. 

More detailed results can be found in HI. 
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